#include "Crossmatch.hh"


int gDepth = 0;

void PrintTabs(int depth){
  for(int i = 0;i< depth;i++) cout<<"  ";
}


//------------------------------------------------------
// * ParseCMFile
//------------------------------------------------------
//
void 
Crossmatch::Read(string FileName)
{
  ifstream infile;
  infile.open(FileName.c_str());
  if (infile.fail()){
	 cout<<"ERROR: Unable to open: "<<FileName<<endl;
	 cout<<"ERROR: Endseq location terminated. "<<endl;
	 exit(1);
  }
  
  string tempLine;
  // This will be the last line before the list of significant hits
  GetNextLineContaining(infile,"Maximal",tempLine);
  
  if (infile.eof()){
	 cout<<"ERROR: Prematurely reached end of file: "<<FileName<<endl;
	 cout<<"ERROR: End seq locator terminated. "<<endl;
	 exit(1);
  }

  // Go through the results lines, extracting important information.   
  bool bEndFound = false;
  int offset = 0,temp;
  CrossmatchResult cmtemp;
  do{
	 GetLine(infile,tempLine);
	 // see if the end of hit list has been reached. 
	 if (!Contains(tempLine,"matching")){
		if (!Contains(tempLine,"Score:")){
		  temp = GetNextInteger(tempLine,(int)0,offset);
		  // check for lines with no numbers (i.e. blank lines between results
		  // if it is a blank line, then do nothing, and just get next line.  
		  if (offset>=0){
			 // Extract all of the interesting bits out of this line
			 cmtemp.mScore = temp;
			 cmtemp.mSubstitutions = (float)GetNextDouble(tempLine,offset,offset);
			 cmtemp.mDeletions = (float)GetNextDouble(tempLine,offset,offset);
			 cmtemp.mInsertions = (float)GetNextDouble(tempLine,offset,offset);
			 cmtemp.mSourceName = GetNextWSDelimitedString(tempLine,offset,offset);
			 cmtemp.mSourceStart = (long)GetNextInteger(tempLine,offset,offset);
			 cmtemp.mSourceEnd = (long)GetNextInteger(tempLine,offset,offset);
			 // number we don't care to save
			 temp = GetNextInteger(tempLine,offset,offset); 
			 offset++; // skip )

			 // Either there will be a C here, for compliment strand, or not.  
			 string sTemp = GetNextWSDelimitedString(tempLine,offset,offset);
			 if (sTemp == "C"){
				cmtemp.mbIsReverse = true;
				cmtemp.mQueryName = GetNextWSDelimitedString(tempLine,offset,offset);
				temp = GetNextInteger(tempLine,offset,offset); // number we don't care to save
				offset++; // skip )
				cmtemp.mQueryStart = GetNextInteger(tempLine,offset,offset);
				cmtemp.mQueryEnd = GetNextInteger(tempLine,offset,offset);
			 }else{
				cmtemp.mbIsReverse = false;
				cmtemp.mQueryName = sTemp;			  
				cmtemp.mQueryStart = GetNextInteger(tempLine,offset,offset);
				cmtemp.mQueryEnd = GetNextInteger(tempLine,offset,offset);
				temp = GetNextInteger(tempLine,offset,offset); // number we don't care to save
				offset++; // skip )
			 }

			 // If there is a * at the end of the line then this hit is partially 
			 // contained in another higher scoring hit. 
			 if (Contains(tempLine,"*")) cmtemp.mbIsContained = true;
			 else cmtemp.mbIsContained = false;

			 push_back(cmtemp);
		  }
		}
	 }else{
		bEndFound = true;
	 }
  }while((!infile.eof()) && (!bEndFound));
}

//----------------------------------------
// * GetSinks
//----------------------------------------
//
void 
Crossmatch::GetSinks(string &SourceName,vector<string> &Sinks,int minScore)
{
  for(int i = 0;i< size();i++){
	 if ((*this)[i].mSourceName == SourceName){
		if ((*this)[i].mScore >= minScore){
		  Sinks.push_back((*this)[i].mQueryName);
		}
	 }
  }
}


//----------------------------------------
// * DFSSink
//----------------------------------------
//
void 
Crossmatch::DFSSink(string &SourceName,
				 StringSet &SomeComponent,
				 StringSet &ThisComponent,
				 int minScore)
{
  gDepth++;

  vector<string> Sinks;
  GetSinks(SourceName,Sinks,minScore);
  
  for(int i = 0;i< Sinks.size();i++){
	 // Don't go if it is in some previous component
	 if (!SomeComponent.Contains(Sinks[i])){
		//PrintTabs(gDepth);
		//cout<<Sinks[i]<<" ADD (DFSSink)"<<endl;
		SomeComponent.Add(Sinks[i]);
		ThisComponent.Add(Sinks[i]);

		DFSSink(Sinks[i],SomeComponent,ThisComponent,minScore);
	 }else{
		//PrintTabs(gDepth);
		//cout<<Sinks[i]<<" in some component (DFSSink)"<<endl;
	 }
  }

  gDepth--;
}

//----------------------------------------
// * AddInverseEdges
//----------------------------------------
// KJD BUG  Right now this is just an approximation. 
// 
// Treated as a graph, this is bipartate (none of sources connect to other sources, 
// none of queries connect to other queries).  Adding the symmetric edges makes
// it possible to do a simple DFS for components. 
void 
Crossmatch::AddInverseEdges()
{
  Crossmatch CMTemp;

  for(int i = 0;i< size();i++){
	 int nextIdx = CMTemp.size();
	 CMTemp.resize(nextIdx+1);
	 CMTemp[nextIdx] = (*this)[i];
	 CMTemp[nextIdx].mSourceName = (*this)[i].mQueryName;
	 CMTemp[nextIdx].mQueryName = (*this)[i].mSourceName;	 
  }

  for(int i = 0;i< CMTemp.size();i++){
	 push_back(CMTemp[i]);
  }
}


//----------------------------------------
// * FindConnectedComponents
//----------------------------------------
//
void 
Crossmatch::FindConnectedComponents(vector<StringSet> &Components,int minScore)
{

  StringSet SomeComponent;
  Crossmatch CMTemp = (*this);

  // This doesn't really add the inverse, just an approximation. 
  CMTemp.AddInverseEdges();  

  for(int i = 0;i< CMTemp.size();i++){
	 if (CMTemp[i].mScore >= minScore){
		string SourceName = CMTemp[i].mSourceName;

		// If this crossmatch entry is not already in some component
		if (!SomeComponent.Contains(SourceName)){
		  //cout<<SourceName<<" ADD (FCC)"<<endl;
		  
		  // Add it to a new connected component
		  int nextComponent = Components.size();
		  Components.resize(nextComponent+1);
		  Components[nextComponent].Add(SourceName);
		  
		  // Record that it is in some component.  
		  SomeComponent.Add(SourceName);
		  
		  // DFS on each of it's children to fill out this connected component. 
		  vector<string> Sinks;
		  CMTemp.GetSinks(SourceName,Sinks,minScore);
		  
		  for(int c = 0; c< Sinks.size();c++){
			 if (!SomeComponent.Contains(Sinks[c])){
				SomeComponent.Add(Sinks[c]);
				Components[nextComponent].Add(Sinks[c]);
				// cout<<Sinks[c]<<"Add child"<<endl;
			 }
			 
			 CMTemp.DFSSink(Sinks[c],SomeComponent,Components[nextComponent],minScore);
		  }
		}else{
		  //		cout<<SourceName<<" in some component (FCC)"<<endl;
		}
	 } // if score >= minScore
  }// For each CMTemp
}


/*

sw-score %substitute  %deletions %insertions  QueryID  QueryStart QueryEnd (#basespastend)  C  SourceID (bases prior) SStart SEnd 


Standard output

  The standard output lists matches between any sequence in the first
input file (the "query" sequences) and a sequence in the second (or
later) input file (the "subject" sequences); or, if a single input
file is provided, the matches between any two sequences in this file.
The matches that are reported are controlled by the command line
options -minscore and -masklevel, as well as by the options that
control scoring of the alignments and the band search (see section IV
above).  The reported matches are ordered by query, and for each query
by the position of the start of the alignment within the query.
  For each reported match, an initial output line gives summary
information:

 Example: 

440  2.38 1.39 0.79  hh44a1.s1       33   536 (    0)  C 00311     ( 3084)  8277   7771  * 

Interpretation:
  440 = smith-waterman score of the match (complexity-adjusted, by default).
  2.38 = %substitutions in matching region
  1.39 = %deletions (in 1st seq rel to 2d) in matching region
  0.79 = %insertions (in 1st seq rel to 2d) in matching region 
  hh44a1.s1 = id of 1st sequence
  33 = starting position of match in 1st sequence
  536 = ending position of match in 1st sequence
  (0) = no. of bases in 1st sequence past the ending position of match
         (so 0 means that the match extended all the way to the end of 
          the 1st sequence)
  C 00311 : match is with the Complement of sequence 00311
  ( 3084) : there are 3084 bases in (complement of) 2d sequence prior to 
        beginning of the match
  8277 = starting position of match in 2d sequence (using top-strand 
         numbering)
  7771 =  ending position of match in 2d sequence
  * indicates that there is a higher-scoring match whose domain partly
includes the domain of this match.

*/
